function fig = PlotMagnetometerSensitivity(obj, sensitivityData)
    if nargin == 1
        sensitivityData = obj.MagnetometorSensitivity;
    end
    
    fList= sensitivityData.PS_DataArray.fList;
    psValue = sensitivityData.PS_DataArray.psValue;
    enbw =  sensitivityData.ENBW;
    coeff =  sensitivityData.Coeff;
    sig_amp = sensitivityData.Signal.amplitude;
    sig_freq = sensitivityData.Signal.frequency;
    
    psd = psValue/sqrt(enbw) / coeff * 1e3;
    ncount = histcounts(psd, 0:1000); % count sensitivity fro 0 ~ 1000 pT/sqrt(Hz)
    [~, loc, wid ] = findpeaks(ncount, 'NPeaks', 1, 'SortStr', 'descend');
    meanSensitivity = mean(psd(psd<loc+3*wid & psd>loc-3*wid));

    fig = figure('Name', sprintf('Magnetometer Sensitivity: [%s]', datestr(sensitivityData.timestamp)), 'WindowStyle', 'Docked'); ax = axes('Parent', fig);
    semilogy(ax, fList, psd, 'b.-', [fList(1), fList(end)], [meanSensitivity, meanSensitivity], 'r--');
    xlabel(ax, 'frequency (Hz)'); ylabel(ax, 'sensitivity pT/sqrt(Hz)');
    xlim(ax, minmax(fList(:)')); 
    title(ax, sprintf('Signal = %3.2f nT @ %2.1f Hz, ENBW = %3.2f Hz', sig_amp, sig_freq, enbw));
    legend(ax, {'PSD', sprintf('Background = %3.2f pT/\\surd{Hz}', meanSensitivity)}, 'Location', 'best', 'Interpreter', 'tex');
    grid on;
end

